Mapping Areas That Could Benefit From Affordable Housing in New Jersey

1.1 Problem Statement

Can a SHINY application be developed that highlights the areas of New Jersey that could most benefit from affordable hsouing? Can users have the ability to select the variable(s) that are used in the screening criteria of the SHINY application?

2.1 Introduction

+INSERT NARRATIVE ON CHALLENGES OF COMMUNITY DEVELOPMENT; STIGMA OF AFFORDABLE HOUSING; LITERATURE REVIEW OF EIGHT VARIABLES BELOW WHERE AFFORDABLE HOUSING CAN BENEFIT A COMMUNITY+

2.2 Literature Review

+INSERT NARRATIVE ON LITERATURE REVIEW+

2.2.1 Affordability

The National Low Income Housing Coalition publishes an annual report that analyzes the levels of housing affordability at the national, state, and county levels. Hourly wage needed to afford a two bedroom apartment - paying no more than 30% of their income towards housing - is one metric that is analyzed. Shortages of affordable rental housing leads to increases in this metric. (Source: National Low Income Housing Coalition)

2.2.2 Property Values

The Center for Housing Policy published a report that noted that affordable housing had a neutral or positive effect on the property values of neighboring parcels. (Source: Center for Housing Policy)

2.2.3 Diversity

The Brookings Institute and The Urban Institute published a report that noted production programs that increase availabilty of affordable housing promotes racial diversity. (Source: The Brookings Institute)

2.2.4 Travel Time

The Brookings Institute publised an article that noted an increase in affordable housing helps reduce commute times as provides lower-income workers increased housing options near places of employment. (Source: The Brookings Institute)

2.2.5 Real Estate Taxes

The Center for Housing Policy published a report that noted that affordable housing can lead to higher real estate income for a municipalit as it often reduces instances of foreclosure which leads to depressed property values and tax collections. (Source: Center for Housing Policy)

2.2.6 Education

The Center for Housing Policy published a report that noted affordable housing may improve educational outcomes of children as it helps increase housing stability reduces family homelessness. (Source: Center for Housing Policy)

2.2.7 Asthma

The American Public Health Assoication published an article that noted the green and sustainable development practices in affordable housing development and renovations can lead to lower asthma rates of residents. (Source: American Public Health Assoication)

2.2.8 Poverty 65+

The Urban Institute published a report that noted affordable housing helps reduce poverty in adults 65 and older as it allows them to age in place which is often a more affordable option than relocating to assisted-living facilities. (Source: The Urban Institute)

2.3 Model and SHINY Application

+INSERT NARRATIVE ON MODEL AND SHINY APPLICATION+

3.1 Methods

A model was developed using R statistical package to aggregate and process all data to create a final csv file that will be used in the SHINY application. Details on data collection, analysis, processing, and aggregation are included below. Code for the SHINY application as well as a URL are also found below.

3.2 Setup

The model begins by downloading the required R libraries. It then sets up plot and map themes as well as a color palette. This allows the model to produce consistent visualizations throughout.

knitr::opts_chunk$set(cache=TRUE)

setwd("C:/Users/scott/Desktop/MUSA800")
options(scipen=999)

library(tidyverse)
library(sf)
library(lubridate)
library(tigris)
library(tidycensus)
library(viridis)
library(riem)
library(gridExtra)
library(knitr)
library(kableExtra)
library(RSocrata)
library(gganimate)
library(ggplot2)
library(dplyr)
library(MASS)
library(QuantPsyc)
library(caret)
library(spatstat)
library(spdep)
library(FNN)
library(grid)
library(RColorBrewer)
library(wesanderson)
library(gridExtra)
library(grid)
library(stats)
library(lme4)
library(rgdal)
library(shiny)

plotTheme <- theme(
  plot.title =element_text(size=12),
  plot.subtitle = element_text(size=8),
  plot.caption = element_text(size = 6),
  axis.text.x = element_text(size = 10, angle = 45, hjust = 1),
  axis.text.y = element_text(size = 10),
  axis.title.y = element_text(size = 10),
  # Set the entire chart region to blank
  panel.background=element_blank(),
  plot.background=element_blank(),
  #panel.border=element_rect(colour="#F0F0F0"),
  # Format the grid
  panel.grid.major=element_line(colour="#D0D0D0",size=.2),
  axis.ticks=element_blank())

mapTheme <- function(base_size = 12) {
  theme(
    text = element_text( color = "black"),
    plot.title = element_text(size = 14,colour = "black"),
    plot.subtitle=element_text(face="italic"),
    plot.caption=element_text(hjust=0),
    axis.ticks = element_blank(),
    panel.background = element_blank(),axis.title = element_blank(),
    axis.text = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=2)
  )
}

palette5 <- c("#eff3ff","#bdd7e7","#6baed6","#3182bd","#08519c")
palette4 <- c("#D2FBD4","#92BCAB","#527D82","#123F5A")
palette2 <- c("#6baed6","#08519c")

3.3 Add Boundaries

State, County, and municipal boundaries are added to the model. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

#-------------------Create NJ Boundary-------------------

US_Boundary <- st_read("C:/Users/scott/Desktop/MUSA800/Data/states.shp") 
US_Boundary <- US_Boundary %>% st_transform(crs=102271)
NJ_Boundary<- US_Boundary[(US_Boundary$STATE_ABBR=="NJ"),]

NJ_Boundary <- 
  NJ_Boundary %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 102271, agr = "constant") 

NJ_Towns <- st_read("https://opendata.arcgis.com/datasets/3d5d1db8a1b34b418c331f4ce1fd0fef_2.geojson")
NJ_Towns <- 
  NJ_Towns %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 102271, agr = "constant") %>%
  st_transform(st_crs(NJ_Boundary))
NJ_Towns <- NJ_Towns %>% dplyr::select("NAME", "MUN_CODE", "geometry")

NJ_Counties <- st_read("https://opendata.arcgis.com/datasets/5f45e1ece6e14ef5866974a7b57d3b95_1.geojson")
NJ_Counties <- 
  NJ_Counties %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 102271, agr = "constant") %>%
  st_transform(st_crs(NJ_Boundary))

ggplot() + 
  geom_sf(data=NJ_Counties, fill=NA, colour="black") +
  mapTheme()

3.4 Create Fishnet

A fishnet was applied to produce equal sized polygons (2,500 sqare meters) over the state and use them as a unit of measure for my analysis. Fishnet and boundaries were converted from crs of 102271 to 4326 so it would work in SHINY application.

knitr::opts_chunk$set(cache=TRUE)

#-------------------Create Fishnet-------------------

createFishnet <- function(x) {
  fishnet <- 
    st_make_grid(NJ_Boundary, cellsize = x) %>%
    st_sf()
  fishnet <- 
    fishnet[NJ_Boundary,] %>%
    mutate(uniqueID = rownames(.)) %>%
    dplyr::select(uniqueID)
  fishnet <- 
    fishnet %>% 
    st_as_sf(coords = c("Longitude", "Latitude"), crs = 102271, agr = "constant") %>%
    st_transform(st_crs(NJ_Boundary))
  fishnet <<- st_intersection(fishnet, NJ_Counties)%>%
    dplyr::select(uniqueID, geometry,COUNTY)}

createFishnet(2500)

#-------------------Convert CRS to 4326 so it will work in SHINY-------------------

US_Boundary <- US_Boundary %>% st_transform(crs=4326)
NJ_Boundary<- US_Boundary[(US_Boundary$STATE_ABBR=="NJ"),]
NJ_Boundary <- 
  NJ_Boundary %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") 

NJ_Towns <- st_read("https://opendata.arcgis.com/datasets/3d5d1db8a1b34b418c331f4ce1fd0fef_2.geojson")
NJ_Towns <- 
  NJ_Towns %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
  st_transform(st_crs(NJ_Boundary))
NJ_Towns <- NJ_Towns %>% dplyr::select("NAME", "MUN_CODE", "geometry")

NJ_Counties <- st_read("https://opendata.arcgis.com/datasets/5f45e1ece6e14ef5866974a7b57d3b95_1.geojson")
NJ_Counties <- 
  NJ_Counties %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
  st_transform(st_crs(NJ_Boundary))

fishnet <- 
  fishnet %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
  st_transform(st_crs(NJ_Boundary))

ggplot() +
  geom_sf(data=NJ_Counties, fill=NA, colour="black") +
  geom_sf(data=fishnet, fill=NA, colour="GRAY") +
  labs(title = "Fishnet in New Jersey") +
  mapTheme()

3.5 Add Data

3.5.1 Affordability Data

The model measures affordability by the hourly wage needed to afford a two bedroom apartment by county. Data was collected from the National Low Income Housing Coalition through their annual Out of Reach reports. County level data was pulled from reports of years 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 and entered into a csv file. This was then inputted into a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the underlying value in each year. Below is a facetted map series of the output.

knitr::opts_chunk$set(cache=TRUE)

#-------------------AFFORDABILITY-------------------
#Hourly wage necessary to afford 2 BR FMR
#https://nlihc.org/sites/default/files/oor/2009-OOR.pdf
#https://nlihc.org/sites/default/files/oor/2010-OOR.pdf
#https://nlihc.org/sites/default/files/oor/2011-OOR.pdf
#https://nlihc.org/sites/default/files/oor/2012-OOR.pdf
#https://reports.nlihc.org/sites/default/files/oor/2013-OOR-NJ_0.pdf
#https://reports.nlihc.org/sites/default/files/oor/2014-OOR-NJ.pdf
#https://nlihc.org/sites/default/files/oor/OOR_2015_FULL.pdf
#https://nlihc.org/sites/default/files/oor/OOR_2016.pdf
#https://reports.nlihc.org/sites/default/files/oor/OOR_2017_0.pdf
#https://nlihc.org/sites/default/files/oor/OOR_2018.pdf


NJ_Afford <- read.csv("C:/Users/scott/Desktop/MUSA800/Data/NJ_Affordability.csv")
names(NJ_Afford)[1] <- "COUNTY"

AFFORD.Net <- fishnet
AFFORD.Net$L2009 <- NJ_Afford$L2009[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)] 
AFFORD.Net$L2010 <- NJ_Afford$L2010[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)] 
AFFORD.Net$L2011 <- NJ_Afford$L2011[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)] 
AFFORD.Net$L2012 <- NJ_Afford$L2012[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)] 
AFFORD.Net$L2013 <- NJ_Afford$L2013[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)] 
AFFORD.Net$L2014 <- NJ_Afford$L2014[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)] 
AFFORD.Net$L2015 <- NJ_Afford$L2015[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)] 
AFFORD.Net$L2016 <- NJ_Afford$L2016[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)] 
AFFORD.Net$L2017 <- NJ_Afford$L2017[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)] 
AFFORD.Net$L2018 <- NJ_Afford$L2018[match(AFFORD.Net$COUNTY, NJ_Afford$COUNTY)] 

##Facett Maps

AFFORD.Map <- AFFORD.Net %>%
  gather(PIS, AFFORD, L2009:L2018)

AFFORD.Map$PIS <- substring(AFFORD.Map$PIS, 2)

ggplot() + 
  geom_sf(data = AFFORD.Map, aes(fill = AFFORD)) + 
  labs(title = "Wage Required-2 Bedroom Apartment, New Jersey") +
  mapTheme() + 
  facet_wrap(~PIS, ncol = 5) + 
  scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")

AFFORD.Map$YEAR.ID <- paste(AFFORD.Map$PIS,AFFORD.Map$uniqueID)

The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Affordability dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which affordability has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Calculate Rate of Change

data <- AFFORD.Net %>%
  dplyr::select(uniqueID,L2009:L2018)

data[is.na(data)] <- 0

data$geometry <- NULL

data$slope <- NA

data[1:nrow(data), "slope"] <- apply(data[1:nrow(data), 
                                          c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
                                     1, function(row_i){lm(unlist(row_i) ~ unlist(data[1, 
                                                                                       c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
                                     })

AFFORD.slope <- fishnet
AFFORD.slope$slope <- data$slope[match(data$uniqueID, AFFORD.slope$uniqueID)]
#RE.slope$slope[is.na(RE.slope$slope)] <- 0

ggplot() + 
  geom_sf(data = AFFORD.slope, aes(fill = slope)) + 
  labs(title = "Slope of Affordability (2009-2018), New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")

A score is then applied to the slope values. This is done by assigning a value of “1” to the highest half of slope values which are those fishnet cells that have the fastest reduction in affordability, and a value of “2” for the lowest half of the slope values which are those fishnet cells that have the slowest reduction in affordability. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145

AFFORD.slope$slope2 <- AFFORD.slope$slope * -1

AFFORD.slope <- within(AFFORD.slope, score <- as.integer(cut(slope2, quantile(slope2, probs=0:2/2),
                                                                  include.lowest=TRUE)))

ggplot() + 
  geom_sf(data = AFFORD.slope, aes(fill = score)) + 
  labs(title = "Affordabilty Score, New Jersey") +
  mapTheme() + 
   scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")

3.5.2 Property Value Data

The model measures property value by the mean property value by census tract. Data was collected from the US Census Bureau using tidycensus for a ten year period (2009-2018) saved to a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the average mean property value in each year. Here the average was taken to account for grid cells that had multiple census tract boudaries. Below is a facetted map series of the output.

knitr::opts_chunk$set(cache=TRUE)

#-------------------PROPERTY VALUE-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference: https://furmancenter.org/files/media/Dont_Put_It_Here.pdf
#https://censusreporter.org/

census_api_key("7c4f7a4ee0b34622dcab1bef185348b87ae7b355", overwrite = TRUE)

NJCensus <- 
  get_acs(geography = "tract", 
          variables = c("B25077_001"),  
          year = 2010, 
          state = "NJ", 
          geometry = TRUE, 
          output = "wide") %>%
  rename(PropVal = B25077_001E) %>%
  dplyr::select(PropVal, GEOID, geometry) 

NJCensus <- 
  NJCensus %>%
  as.data.frame() %>%
  distinct(GEOID, .keep_all = TRUE) %>%
  #dplyr::select(GEOID, geometry) %>% 
  st_sf

NJCensus <- 
  NJCensus %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
  st_transform(st_crs(NJ_Boundary), 
               join=st_intersects,
               left = TRUE) 

census <- function(x) {
  NJCensus <<- 
    get_acs(geography = "tract", 
            variables = c("B25077_001"),  
            year = x, 
            state = "NJ", 
            geometry = TRUE, 
            output = "wide") %>%
    rename(PropVal = B25077_001E) %>%
    dplyr::select(PropVal, GEOID, geometry) %>%
    as.data.frame() %>%
    distinct(GEOID, .keep_all = TRUE) %>%
    #dplyr::select(GEOID, geometry) %>% 
    st_sf%>%
    st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(NJ_Boundary), 
                 join=st_intersects,
                 left = TRUE) 
  NJCensus$PIS <<-x
  
}

AVG <- function(x) {
  NET <- NET[NET$PIS <= x,]
  AVG <- aggregate( PropVal ~ uniqueID, NET, mean)  
  PROPVAL.Net$y <<- AVG$PropVal[match(PROPVAL.Net$uniqueID, AVG$uniqueID)] 
  PROPVAL.Net$y[is.na(PROPVAL.Net$y)] <<- 0
}

PROPVAL.Net <- fishnet
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2009"
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2010"
census(2011); NET <- st_intersection(fishnet, NJCensus)
AVG(2011); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2011"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2012"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2013"
census(2014); NET <- st_intersection(fishnet, NJCensus)
AVG(2014); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2014"
census(2015); NET <- st_intersection(fishnet, NJCensus)
AVG(2015); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2015"
census(2016); NET <- st_intersection(fishnet, NJCensus)
AVG(2016); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2016"
census(2017); NET <- st_intersection(fishnet, NJCensus)
AVG(2017); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2017"
census(2018); NET <- st_intersection(fishnet, NJCensus)
AVG(2018); colnames(PROPVAL.Net)[colnames(PROPVAL.Net)=="y"] <- "L2018"


###Facett Maps

PROPVAL.Map <- PROPVAL.Net %>%
  gather(PIS, PROPVAL, L2009:L2018)

PROPVAL.Map$PIS <- substring(PROPVAL.Map$PIS, 2)

ggplot() + 
  geom_sf(data = PROPVAL.Map, aes(fill = PROPVAL)) + 
  labs(title = "Property Value, New Jersey") +
  mapTheme() + 
  facet_wrap(~PIS, ncol = 5) + 
  scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")

PROPVAL.Map$YEAR.ID <- paste(PROPVAL.Map$PIS,PROPVAL.Map$uniqueID)

The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Property Value dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which mean property value has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Calculate Rate of Change

data <- PROPVAL.Net %>%
  dplyr::select(uniqueID,L2009:L2018)

data$geometry <- NULL

data$slope <- NA

data[1:nrow(data), "slope"] <- apply(data[1:nrow(data), 
                                          c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
                                     1, function(row_i){lm(unlist(row_i) ~ unlist(data[1, 
                                                                                       c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
                                     })

PROPVAL.slope <- fishnet
PROPVAL.slope$slope <- data$slope[match(data$uniqueID, PROPVAL.slope$uniqueID)]
#PROPVAL.slope$slope[is.na(PROPVAL.slope$slope)] <- 0

ggplot() + 
  geom_sf(data = PROPVAL.slope, aes(fill = slope)) + 
  labs(title = "Slope of Property Value (2009-2018), New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")

A score is then applied to the slope values. This is done by assigning a value of “1” to the lowest half of slope values which are those fishnet cells that have the slowest increase in property values, and a value of “2” for the highest half of the slope values which are those fishnet cells that have the fastest increase in property values. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145

#PROPVAL.slope$slope2 <- PROPVAL.slope$slope * -1

PROPVAL.slope <- within(PROPVAL.slope, score <- as.integer(cut(slope, quantile(slope, probs=0:2/2),
                                                             include.lowest=TRUE)))

ggplot() + 
  geom_sf(data = PROPVAL.slope, aes(fill = score)) + 
  labs(title = "Property Value Score, New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")

3.5.3 Diversity Data

The model measures diversity by the percent white of population by census tract. Data was collected from the US Census Bureau using tidycensus for a ten year period (2009-2018) saved to a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the average percent white in each year. Here the average was taken to account for grid cells that had multiple census tract boudaries. Below is a facetted map series of the output.

knitr::opts_chunk$set(cache=TRUE)

#-------------------PERCENT WHITE-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference: https://www.brookings.edu/wp-content/uploads/2016/06/housingreview.pdf
#https://censusreporter.org/

NJCensus <- 
  get_acs(geography = "tract", 
          variables = c("B01003_001","B02001_002"),  
          year = 2011, 
          state = "NJ", 
          geometry = TRUE, 
          output = "wide") %>%
  rename(Total_Pop =  B01003_001E, White_Pop = B02001_002E) %>%
  dplyr::select(Total_Pop, White_Pop, GEOID, geometry) %>%
  mutate(Percent_White = White_Pop / Total_Pop)

NJCensus <- 
  NJCensus %>%
  as.data.frame() %>%
  distinct(GEOID, .keep_all = TRUE) %>%
  #dplyr::select(GEOID, geometry) %>% 
  st_sf

NJCensus <- 
  NJCensus %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
  st_transform(st_crs(NJ_Boundary), 
               join=st_intersects,
               left = TRUE) 

census <- function(x) {
  NJCensus <<- 
    get_acs(geography = "tract", 
            variables = c("B01003_001","B02001_002"),  
            year = x, 
            state = "NJ", 
            geometry = TRUE, 
            output = "wide") %>%
    rename(Total_Pop =  B01003_001E, White_Pop = B02001_002E) %>%
    dplyr::select(Total_Pop, White_Pop, GEOID, geometry) %>%
    mutate(Percent_White = White_Pop / Total_Pop)%>%
    as.data.frame() %>%
    distinct(GEOID, .keep_all = TRUE) %>%
    #dplyr::select(GEOID, geometry) %>% 
    st_sf%>%
    st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(NJ_Boundary), 
                 join=st_intersects,
                 left = TRUE) 
  NJCensus$PIS <<-x
  
}

AVG <- function(x) {
  NET <- NET[NET$PIS <= x,]
  AVG <- aggregate( Percent_White ~ uniqueID, NET, mean)  
  PCTWHITE.Net$y <<- AVG$Percent_White[match(PCTWHITE.Net$uniqueID, AVG$uniqueID)] 
  PCTWHITE.Net$y[is.na(PCTWHITE.Net$y)] <<- 0
}

PCTWHITE.Net <- fishnet
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2009"
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2010"
census(2011); NET <- st_intersection(fishnet, NJCensus)
AVG(2011); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2011"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2012"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2013"
census(2014); NET <- st_intersection(fishnet, NJCensus)
AVG(2014); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2014"
census(2015); NET <- st_intersection(fishnet, NJCensus)
AVG(2015); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2015"
census(2016); NET <- st_intersection(fishnet, NJCensus)
AVG(2016); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2016"
census(2017); NET <- st_intersection(fishnet, NJCensus)
AVG(2017); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2017"
census(2018); NET <- st_intersection(fishnet, NJCensus)
AVG(2018); colnames(PCTWHITE.Net)[colnames(PCTWHITE.Net)=="y"] <- "L2018"


##Facett Maps

PCTWHITE.Map <- PCTWHITE.Net %>%
  gather(PIS, PCTWHITE, L2009:L2018)

PCTWHITE.Map$PIS <- substring(PCTWHITE.Map$PIS, 2)

ggplot() + 
  geom_sf(data = PCTWHITE.Map, aes(fill = PCTWHITE)) + 
  labs(title = "Percent White, New Jersey") +
  mapTheme() + 
  facet_wrap(~PIS, ncol = 5) + 
  scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")

PCTWHITE.Map$YEAR.ID <- paste(PCTWHITE.Map$PIS,PCTWHITE.Map$uniqueID)

The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Diversity dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which percent white has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Calculate Rate of Change

data <- PCTWHITE.Net %>%
  dplyr::select(uniqueID,L2009:L2018)

data$geometry <- NULL

data$slope <- NA

data[1:nrow(data), "slope"] <- apply(data[1:nrow(data), 
                                          c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
                                     1, function(row_i){lm(unlist(row_i) ~ unlist(data[1, 
                                                                                       c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
                                     })

PCTWHITE.slope <- fishnet
PCTWHITE.slope$slope <- data$slope[match(data$uniqueID, PCTWHITE.slope$uniqueID)]
#PCTWHITE.slope$slope[is.na(PCTWHITE.slope$slope)] <- 0

ggplot() + 
  geom_sf(data = PCTWHITE.slope, aes(fill = slope)) + 
  labs(title = "Slope of Diversity (2009-2018), New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")

A score is then applied to the slope values. This is done by assigning a value of “1” to the highest half of slope values which are those fishnet cells that have the fastest increase in percent white, and a value of “2” for the lowest half of the slope values which are those fishnet cells that have the fastest decrease in percent white. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145

PCTWHITE.slope$slope2 <- PCTWHITE.slope$slope * -1

PCTWHITE.slope <- within(PCTWHITE.slope, score <- as.integer(cut(slope2, quantile(slope2, probs=0:2/2),
                                                               include.lowest=TRUE)))

ggplot() + 
  geom_sf(data = PCTWHITE.slope, aes(fill = score)) + 
  labs(title = "Diversity Score, New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")

3.5.4 Travel Time Data

The model measures travel time by percent of population that has over 90 minute commute to work by census tract. Data was collected from the US Census Bureau using tidycensus for a ten year period (2009-2018) saved to a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the average travel time in each year. Here the average was taken to account for grid cells that had multiple census tract boudaries. Below is a facetted map series of the output.

knitr::opts_chunk$set(cache=TRUE)

#-------------------TRAVEL TIME-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference:  http://www.reconnectingamerica.org/assets/Uploads/20120904AHpolicybrief.pdf
#https://censusreporter.org/

NJCensus <- 
  get_acs(geography = "tract", 
          variables = c("B01003_001","B08303_013"),  
          year = 2009, 
          state = "NJ", 
          geometry = TRUE, 
          output = "wide") %>%
  rename(Total_Pop =  B01003_001E, TRAVEL = B08303_013E) %>%
  dplyr::select(Total_Pop, TRAVEL, GEOID, geometry) %>%
  mutate(AVGCOMMUTE = TRAVEL / Total_Pop)

NJCensus <- 
  NJCensus %>%
  as.data.frame() %>%
  distinct(GEOID, .keep_all = TRUE) %>%
  #dplyr::select(GEOID, geometry) %>% 
  st_sf

NJCensus <- 
  NJCensus %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
  st_transform(st_crs(NJ_Boundary), 
               join=st_intersects,
               left = TRUE) 


census <- function(x) {
  NJCensus <<- 
    get_acs(geography = "tract", 
            variables = c("B01003_001","B08303_013"),  
            year = x, 
            state = "NJ", 
            geometry = TRUE, 
            output = "wide") %>%
    rename(Total_Pop =  B01003_001E, TRAVEL = B08303_013E) %>%
    dplyr::select(Total_Pop, TRAVEL, GEOID, geometry) %>%
    mutate(AVGCOMMUTE = TRAVEL / Total_Pop) %>%
    as.data.frame() %>%
    distinct(GEOID, .keep_all = TRUE) %>%
    #dplyr::select(GEOID, geometry) %>% 
    st_sf%>%
    st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(NJ_Boundary), 
                 join=st_intersects,
                 left = TRUE) 
  NJCensus$PIS <<-x
  
  
}

AVG <- function(x) {
  NET <- NET[NET$PIS <= x,]
  AVG <- aggregate( AVGCOMMUTE ~ uniqueID, NET, mean)  
  TRAVEL.Net$y <<- AVG$AVGCOMMUTE[match(TRAVEL.Net$uniqueID, AVG$uniqueID)] 
  TRAVEL.Net$y[is.na(TRAVEL.Net$y)] <<- 0
}

TRAVEL.Net <- fishnet
census(2009); NET <- st_intersection(fishnet, NJCensus)
AVG(2009); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2009"
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2010"
census(2011); NET <- st_intersection(fishnet, NJCensus)
AVG(2011); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2011"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2012"
census(2013); NET <- st_intersection(fishnet, NJCensus)
AVG(2013); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2013"
census(2014); NET <- st_intersection(fishnet, NJCensus)
AVG(2014); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2014"
census(2015); NET <- st_intersection(fishnet, NJCensus)
AVG(2015); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2015"
census(2016); NET <- st_intersection(fishnet, NJCensus)
AVG(2016); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2016"
census(2017); NET <- st_intersection(fishnet, NJCensus)
AVG(2017); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2017"
census(2018); NET <- st_intersection(fishnet, NJCensus)
AVG(2018); colnames(TRAVEL.Net)[colnames(TRAVEL.Net)=="y"] <- "L2018"

##Facett Maps

TRAVEL.Map <- TRAVEL.Net %>%
  gather(PIS, TRAVEL, L2009:L2018)

TRAVEL.Map$PIS <- substring(TRAVEL.Map$PIS, 2)

ggplot() + 
  geom_sf(data = TRAVEL.Map, aes(fill = TRAVEL)) + 
  labs(title = "Travel Time, New Jersey") +
  mapTheme() + 
  facet_wrap(~PIS, ncol = 5) + 
  scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")

TRAVEL.Map$YEAR.ID <- paste(TRAVEL.Map$PIS,TRAVEL.Map$uniqueID)

The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Travel Time dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which travel time has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Calculate Rate of Change

data <- TRAVEL.Net %>%
  dplyr::select(uniqueID,L2009:L2018)

data$geometry <- NULL

data$slope <- NA

data[1:nrow(data), "slope"] <- apply(data[1:nrow(data), 
  c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
  1, function(row_i){lm(unlist(row_i) ~ unlist(data[1, 
  c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
  })

TRAVEL.slope <- fishnet
TRAVEL.slope$slope <- data$slope[match(data$uniqueID, TRAVEL.slope$uniqueID)]
#TRAVEL.slope$slope[is.na(TRAVEL.slope$slope)] <- 0

ggplot() + 
  geom_sf(data = TRAVEL.slope, aes(fill = slope)) + 
  labs(title = "Slope of Travel Time (2009-2018), New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")

A score is then applied to the slope values. This is done by assigning a value of “1” to the highest half of slope values which are those fishnet cells that have the fastest increase in travel, and a value of “2” for the lowest half of the slope values which are those fishnet cells that have the fastest decrease in travel. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145

#TRAVEL.slope$slope2 <- TRAVEL.slope$slope * -1

TRAVEL.slope <- within(TRAVEL.slope, score <- as.integer(cut(slope, quantile(slope, probs=0:2/2),
                                                                 include.lowest=TRUE)))

ggplot() + 
  geom_sf(data = TRAVEL.slope, aes(fill = score)) + 
  labs(title = "Travel Score, New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")

3.5.5 Real Estate Tax Data

The model measures real estate tax by the annual real estate tax collection by municipality. Data was collected from the NJ Department of Community Affairs for a ten year period (2009-2018) and then saved to a csv file. This was then inputted into a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the average mean real estate tax collections in each year. Here the average was taken to account for grid cells that had multiple municipality boudaries. Below is a facetted map series of the output.

knitr::opts_chunk$set(cache=TRUE)

#-------------------RE TAXES-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference: https://providencehousing.org/wp-content/uploads/2014/03/Housing-and-Economic-Development-Report-2011.pdf
#https://censusreporter.org/


#Reference: https://providencehousing.org/wp-content/uploads/2014/03/Housing-and-Economic-Development-Report-2011.pdf
##Data Source: https://www.state.nj.us/dca/divisions/dlgs/resources/property_tax.html
NJ_RETaxes <- read.csv("C:/Users/scott/Desktop/MUSA800/Data/RE_Taxes.csv")

NJ_RETaxes <- NJ_RETaxes %>% rename(NAME = Town)
NJ_RETaxes <- merge(NJ_RETaxes,NJ_Towns,by="NAME")
NJ_RETaxes$X <- NULL
NJ_RETaxes$X.1 <- NULL
NJ_RETaxes[1] <- NULL

NJ_RETaxes <- 
  NJ_RETaxes %>% 
  st_as_sf(crs = 4326, agr = "constant") %>%
  st_transform(st_crs(NJ_Boundary))

NJ_RETaxes.Net <- st_intersection(fishnet, NJ_RETaxes)

RE.Net <- fishnet
AggTax <- aggregate( L2009 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2009 <- AggTax$L2009[match(RE.Net$uniqueID, AggTax$uniqueID)] 
AggTax <- aggregate( L2010 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2010 <- AggTax$L2010[match(RE.Net$uniqueID, AggTax$uniqueID)] 
AggTax <- aggregate( L2011 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2011 <- AggTax$L2011[match(RE.Net$uniqueID, AggTax$uniqueID)] 
AggTax <- aggregate( L2012 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2012 <- AggTax$L2012[match(RE.Net$uniqueID, AggTax$uniqueID)] 
AggTax <- aggregate( L2013 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2013 <- AggTax$L2013[match(RE.Net$uniqueID, AggTax$uniqueID)] 
AggTax <- aggregate( L2014 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2014 <- AggTax$L2014[match(RE.Net$uniqueID, AggTax$uniqueID)] 
AggTax <- aggregate( L2015 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2015 <- AggTax$L2015[match(RE.Net$uniqueID, AggTax$uniqueID)] 
AggTax <- aggregate( L2016 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2016 <- AggTax$L2016[match(RE.Net$uniqueID, AggTax$uniqueID)] 
AggTax <- aggregate( L2017 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2017 <- AggTax$L2017[match(RE.Net$uniqueID, AggTax$uniqueID)] 
AggTax <- aggregate( L2018 ~ uniqueID, NJ_RETaxes.Net, mean)
RE.Net$L2018 <- AggTax$L2018[match(RE.Net$uniqueID, AggTax$uniqueID)] 

##Facett Maps

RE.Map <- RE.Net %>%
  gather(PIS, TAXES, L2009:L2018)

RE.Map$PIS <- substring(RE.Map$PIS, 2)

ggplot() + 
  geom_sf(data = RE.Map, aes(fill = TAXES)) + 
  labs(title = "Real Estate Taxes, New Jersey") +
  mapTheme() + 
  facet_wrap(~PIS, ncol = 5) + 
  scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")

RE.Map$YEAR.ID <- paste(RE.Map$PIS,RE.Map$uniqueID)

The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Real Estate Tax dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which real estate tax collection has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Calculate Rate of Change

data <- RE.Net %>%
  dplyr::select(uniqueID,L2009:L2018)

data[is.na(data)] <- 0

data$geometry <- NULL

data$slope <- NA

data[1:nrow(data), "slope"] <- apply(data[1:nrow(data), 
                                          c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
                                     1, function(row_i){lm(unlist(row_i) ~ unlist(data[1, 
                                                                                       c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
                                     })

RE.slope <- fishnet
RE.slope$slope <- data$slope[match(data$uniqueID, RE.slope$uniqueID)]
#RE.slope$slope[is.na(RE.slope$slope)] <- 0

ggplot() + 
  geom_sf(data = RE.slope, aes(fill = slope)) + 
  labs(title = "Slope of Real Estate Taxes (2009-2018), New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")

A score is then applied to the slope values. This is done by assigning a value of “1” to the lowest half of slope values which are those fishnet cells that have the slowest increase in real estate tax collections, and a value of “2” for the highest half of the slope values which are those fishnet cells that have the fastest increase in real estate tax collections. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145

#RE.slope$slope2 <- RE.slope$slope * -1

RE.slope <- within(RE.slope, score <- as.integer(cut(slope, quantile(slope, probs=0:2/2),
                                                                   include.lowest=TRUE)))

ggplot() + 
  geom_sf(data = RE.slope, aes(fill = score)) + 
  labs(title = "Real Estate Tax Score, New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")

3.5.6 Education Data

The model measures education by high school graduation by census tract. Data was collected from the US Census Bureau using tidycensus for a ten year period (2009-2018) saved to a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the average high school graduation in each year. Here the average was taken to account for grid cells that had multiple census tract boudaries. Below is a facetted map series of the output.

knitr::opts_chunk$set(cache=TRUE)

#-------------------EDUCATION-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference: https://www.nhc.org/wp-content/uploads/2017/03/The-Impacts-of-Affordable-Housing-on-Education-1.pdf
#https://censusreporter.org/


NJCensus <- 
  get_acs(geography = "tract", 
          variables = c("B23006_009"),  
          year = 2009, 
          state = "NJ", 
          geometry = TRUE, 
          output = "wide") %>%
  rename(HIGHSCHOOL = B23006_009E) %>%
  dplyr::select(HIGHSCHOOL, GEOID, geometry) 

NJCensus <- 
  NJCensus %>%
  as.data.frame() %>%
  distinct(GEOID, .keep_all = TRUE) %>%
  #dplyr::select(GEOID, geometry) %>% 
  st_sf

NJCensus <- 
  NJCensus %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
  st_transform(st_crs(NJ_Boundary), 
               join=st_intersects,
               left = TRUE) 

census <- function(x) {
  NJCensus <<- 
    get_acs(geography = "tract", 
            variables = c("B23006_009"),  
            year = x, 
            state = "NJ", 
            geometry = TRUE, 
            output = "wide") %>%
    rename(HIGHSCHOOL = B23006_009E) %>%
    dplyr::select(HIGHSCHOOL, GEOID, geometry) %>%
    as.data.frame() %>%
    distinct(GEOID, .keep_all = TRUE) %>%
    #dplyr::select(GEOID, geometry) %>% 
    st_sf%>%
    st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(NJ_Boundary), 
                 join=st_intersects,
                 left = TRUE) 
  NJCensus$PIS <<-x
  
}

AVG <- function(x) {
  NET <- NET[NET$PIS <= x,]
  AVG <- aggregate( HIGHSCHOOL ~ uniqueID, NET, mean)  
  HIGHSCHOOL.Net$y <<- AVG$HIGHSCHOOL[match(HIGHSCHOOL.Net$uniqueID, AVG$uniqueID)] 
  HIGHSCHOOL.Net$y[is.na(HIGHSCHOOL.Net$y)] <<- 0
}

HIGHSCHOOL.Net <- fishnet
census(2009); NET <- st_intersection(fishnet, NJCensus)
AVG(2009); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2009"
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2010"
census(2011); NET <- st_intersection(fishnet, NJCensus)
AVG(2011); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2011"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2012"
census(2013); NET <- st_intersection(fishnet, NJCensus)
AVG(2013); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2013"
census(2014); NET <- st_intersection(fishnet, NJCensus)
AVG(2014); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2014"
census(2015); NET <- st_intersection(fishnet, NJCensus)
AVG(2015); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2015"
census(2016); NET <- st_intersection(fishnet, NJCensus)
AVG(2016); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2016"
census(2017); NET <- st_intersection(fishnet, NJCensus)
AVG(2017); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2017"
census(2018); NET <- st_intersection(fishnet, NJCensus)
AVG(2018); colnames(HIGHSCHOOL.Net)[colnames(HIGHSCHOOL.Net)=="y"] <- "L2018"


###Facett Maps

HIGHSCHOOL.Map <- HIGHSCHOOL.Net %>%
  gather(PIS, HIGHSCHOOL, L2009:L2018)

HIGHSCHOOL.Map$PIS <- substring(HIGHSCHOOL.Map$PIS, 2)

ggplot() + 
  geom_sf(data = HIGHSCHOOL.Map, aes(fill = HIGHSCHOOL)) + 
  labs(title = "High School Degree, New Jersey") +
  mapTheme() + 
  facet_wrap(~PIS, ncol = 5) + 
  scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")

HIGHSCHOOL.Map$YEAR.ID <- paste(HIGHSCHOOL.Map$PIS,HIGHSCHOOL.Map$uniqueID)

The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Education dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which high school graduation has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Calculate Rate of Change

data <- HIGHSCHOOL.Net %>%
  dplyr::select(uniqueID,L2009:L2018)

data$geometry <- NULL

data$slope <- NA

data[1:nrow(data), "slope"] <- apply(data[1:nrow(data), 
                                          c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
                                     1, function(row_i){lm(unlist(row_i) ~ unlist(data[1, 
                                                                                       c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
                                     })

HIGHSCHOOL.slope <- fishnet
HIGHSCHOOL.slope$slope <- data$slope[match(data$uniqueID, HIGHSCHOOL.slope$uniqueID)]
#PROPVAL.slope$slope[is.na(PROPVAL.slope$slope)] <- 0

ggplot() + 
  geom_sf(data = HIGHSCHOOL.slope, aes(fill = slope)) + 
  labs(title = "Slope of Education (2009-2018), New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")

A score is then applied to the slope values. This is done by assigning a value of “1” to the lowest half of slope values which are those fishnet cells that have the slowest increase in high school graduation, and a value of “2” for the highest half of the slope values which are those fishnet cells that have the fastest increase in high school graduation. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145

#HIGHSCHOOL.slope$slope2 <- HIGHSCHOOL.slope$slope * -1

HIGHSCHOOL.slope <- within(HIGHSCHOOL.slope, score <- as.integer(cut(slope, quantile(slope, probs=0:2/2),
                                                             include.lowest=TRUE)))

ggplot() + 
  geom_sf(data = HIGHSCHOOL.slope, aes(fill = score)) + 
  labs(title = "Education Score, New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")

3.5.7 Asthma Data

The model Asthma by Asthma hospitalizations by County. Data was collected from the NJ Department of Health for a ten year period (2009-2018) and saved to a csv file. This data was then inputted into a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the average travel time in each year. Here the average was taken to account for grid cells that had multiple census tract boudaries. Below is a facetted map series of the output.

knitr::opts_chunk$set(cache=TRUE)

#-------------------ASTHMA-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3000722/

#Source:https://www.nj.gov/health/fhs/chronic/documents/asthma_profiles/asthma_overview.pdf
#Source: https://www-doh.state.nj.us/doh-shad/query/selection/ub/UBSelection.html
#Age-adjusted Rates (Discharges Per 10,000 Standard Population) for Counties-AMI -- Primary Diagnosis

NJ_Asthma <- read.csv("C:/Users/scott/Desktop/MUSA800/Data/NJ_Asthma.csv")
names(NJ_Asthma)[1] <- "COUNTY"

ASTHMA.Net <- fishnet
ASTHMA.Net$L2009 <- NJ_Asthma$L2009[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)] 
ASTHMA.Net$L2010 <- NJ_Asthma$L2010[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)] 
ASTHMA.Net$L2011 <- NJ_Asthma$L2011[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)] 
ASTHMA.Net$L2012 <- NJ_Asthma$L2012[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)] 
ASTHMA.Net$L2013 <- NJ_Asthma$L2013[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)] 
ASTHMA.Net$L2014 <- NJ_Asthma$L2014[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)] 
ASTHMA.Net$L2015 <- NJ_Asthma$L2015[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)] 
ASTHMA.Net$L2016 <- NJ_Asthma$L2016[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)] 
ASTHMA.Net$L2017 <- NJ_Asthma$L2017[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)] 
ASTHMA.Net$L2018 <- NJ_Asthma$L2018[match(ASTHMA.Net$COUNTY, NJ_Asthma$COUNTY)] 

##Facett Maps

ASTHMA.Map <- ASTHMA.Net %>%
  gather(PIS, ASTHMA, L2009:L2018)

ASTHMA.Map$PIS <- substring(ASTHMA.Map$PIS, 2)

ggplot() + 
  geom_sf(data = ASTHMA.Map, aes(fill = ASTHMA)) + 
  labs(title = "Asthma Hospitalizations, New Jersey") +
  mapTheme() + 
  facet_wrap(~PIS, ncol = 5) + 
  scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")

ASTHMA.Map$YEAR.ID <- paste(ASTHMA.Map$PIS,ASTHMA.Map$uniqueID)

The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Asthma dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which asthma hospitalizations has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Calculate Rate of Change

data <- ASTHMA.Net %>%
  dplyr::select(uniqueID,L2009:L2018)

data[is.na(data)] <- 0

data$geometry <- NULL

data$slope <- NA

data[1:nrow(data), "slope"] <- apply(data[1:nrow(data), 
                                          c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
                                     1, function(row_i){lm(unlist(row_i) ~ unlist(data[1, 
                                                                                       c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
                                     })

ASTHMA.slope <- fishnet
ASTHMA.slope$slope <- data$slope[match(data$uniqueID, RE.slope$uniqueID)]
#RE.slope$slope[is.na(RE.slope$slope)] <- 0

ggplot() + 
  geom_sf(data = ASTHMA.slope, aes(fill = slope)) + 
  labs(title = "Slope of Asthma Hospitalizations (2009-2018), New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")

A score is then applied to the slope values. This is done by assigning a value of “1” to the highest half of slope values which are those fishnet cells that have the fastest increase in asthma hospitalizations, and a value of “2” for the lowest half of the slope values which are those fishnet cells that have the fastest decrease in asthma hospitalizations. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145

#ASTHMA.slope$slope2 <- ASTHMA.slope$slope * -1

ASTHMA.slope <- within(ASTHMA.slope, score <- as.integer(cut(slope, quantile(slope, probs=0:2/2),
                                                     include.lowest=TRUE)))

ggplot() + 
  geom_sf(data = ASTHMA.slope, aes(fill = score)) + 
  labs(title = "Asthma Score, New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")

3.5.8 Poverty 65+ Data

The model measures poverty 65+ by poverty status of households 65 and older by census tract. Data was collected from the US Census Bureau using tidycensus for a ten year period (2009-2018) saved to a dataframe with a seperate column for each year. The data was then merged with the fishnet so that every grid cell took on the poverty 65+ in each year. Here the average was taken to account for grid cells that had multiple census tract boudaries. Below is a facetted map series of the output.

knitr::opts_chunk$set(cache=TRUE)

#-------------------SENIORS-------------------
#Enterprise Report: https://homeforallsmc.org/wp-content/uploads/2017/05/Impact-of-Affordable-Housing-on-Families-and-Communities.pdf
#Reference: https://www.urban.org/research/publication/housing-platform-improving-outcomes-older-renters
#https://censusreporter.org/


NJCensus <- 
  get_acs(geography = "tract", 
          variables = c("B17017_008"),  
          year = 2009, 
          state = "NJ", 
          geometry = TRUE, 
          output = "wide") %>%
  rename(POVERTY65 = B17017_008E) %>%
  dplyr::select(POVERTY65, GEOID, geometry) 

NJCensus <- 
  NJCensus %>%
  as.data.frame() %>%
  distinct(GEOID, .keep_all = TRUE) %>%
  #dplyr::select(GEOID, geometry) %>% 
  st_sf

NJCensus <- 
  NJCensus %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
  st_transform(st_crs(NJ_Boundary), 
               join=st_intersects,
               left = TRUE) 

census <- function(x) {
  NJCensus <<- 
    get_acs(geography = "tract", 
            variables = c("B17017_008"),  
            year = x, 
            state = "NJ", 
            geometry = TRUE, 
            output = "wide") %>%
    rename(POVERTY65 = B17017_008E) %>%
    dplyr::select(POVERTY65, GEOID, geometry) %>%
    as.data.frame() %>%
    distinct(GEOID, .keep_all = TRUE) %>%
    #dplyr::select(GEOID, geometry) %>% 
    st_sf%>%
    st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(NJ_Boundary), 
                 join=st_intersects,
                 left = TRUE) 
  NJCensus$PIS <<-x
  
}

AVG <- function(x) {
  NET <- NET[NET$PIS <= x,]
  AVG <- aggregate( POVERTY65 ~ uniqueID, NET, mean)  
  POVERTY65.Net$y <<- AVG$POVERTY65[match(POVERTY65.Net$uniqueID, AVG$uniqueID)] 
  POVERTY65.Net$y[is.na(POVERTY65.Net$y)] <<- 0
}

POVERTY65.Net <- fishnet
census(2009); NET <- st_intersection(fishnet, NJCensus)
AVG(2009); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2009"
census(2010); NET <- st_intersection(fishnet, NJCensus)
AVG(2010); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2010"
census(2011); NET <- st_intersection(fishnet, NJCensus)
AVG(2011); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2011"
census(2012); NET <- st_intersection(fishnet, NJCensus)
AVG(2012); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2012"
census(2013); NET <- st_intersection(fishnet, NJCensus)
AVG(2013); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2013"
census(2014); NET <- st_intersection(fishnet, NJCensus)
AVG(2014); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2014"
census(2015); NET <- st_intersection(fishnet, NJCensus)
AVG(2015); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2015"
census(2016); NET <- st_intersection(fishnet, NJCensus)
AVG(2016); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2016"
census(2017); NET <- st_intersection(fishnet, NJCensus)
AVG(2017); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2017"
census(2018); NET <- st_intersection(fishnet, NJCensus)
AVG(2018); colnames(POVERTY65.Net)[colnames(POVERTY65.Net)=="y"] <- "L2018"


###Facett Maps

POVERTY65.Map <- POVERTY65.Net %>%
  gather(PIS, POVERTY65, L2009:L2018)

POVERTY65.Map$PIS <- substring(POVERTY65.Map$PIS, 2)

ggplot() + 
  geom_sf(data = POVERTY65.Map, aes(fill = POVERTY65)) + 
  labs(title = "Poverty Households 65+, New Jersey") +
  mapTheme() + 
  facet_wrap(~PIS, ncol = 5) + 
  scale_fill_gradient2("Value", low = "#762A83", mid = "white", high = "#1B7837")

POVERTY65.Map$YEAR.ID <- paste(POVERTY65.Map$PIS,POVERTY65.Map$uniqueID)

The model then calculates the rate of change for each fishnet grid cell. This is done by looping an lm() function through the Poverty 65+ dataset which finds the slope of values for years 2009-2018 for each unique ID. The output shows the rate at which poverty status 65+ has changed over time and which direction (increase/decrease) that change is moving. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Calculate Rate of Change

data <- POVERTY65.Net %>%
  dplyr::select(uniqueID,L2009:L2018)

data$geometry <- NULL

data$slope <- NA

data[1:nrow(data), "slope"] <- apply(data[1:nrow(data), 
                                          c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")],
                                     1, function(row_i){lm(unlist(row_i) ~ unlist(data[1, 
                                                                                       c("L2009", "L2010", "L2011", "L2012", "L2013", "L2014", "L2015", "L2016", "L2017", "L2018")]))$coefficients[2]
                                     })

POVERTY65.slope <- fishnet
POVERTY65.slope$slope <- data$slope[match(data$uniqueID, POVERTY65.slope$uniqueID)]
#PROPVAL.slope$slope[is.na(PROPVAL.slope$slope)] <- 0

ggplot() + 
  geom_sf(data = POVERTY65.slope, aes(fill = slope)) + 
  labs(title = "Slope of Poverty Households 65+ (2009-2018), New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Slope", low = "#762A83", mid = "white", high = "#1B7837")

A score is then applied to the slope values. This is done by assigning a value of “1” to the highest half of slope values which are those fishnet cells that have the fastest increase in poverty status 65+, and a value of “2” for the lowest half of the slope values which are those fishnet cells that have the fastest decrease in poverty status 65+. Below is a map of the output.

knitr::opts_chunk$set(cache=TRUE)

###Assign Score
#https://stackoverflow.com/questions/7508229/how-to-create-a-column-with-a-quartile-rank/33896145

#POVERTY65.slope$slope2 <- POVERTY65.slope$slope * -1

POVERTY65.slope <- within(POVERTY65.slope, score <- as.integer(cut(slope, quantile(slope, probs=0:2/2),
                                                                     include.lowest=TRUE)))

ggplot() + 
  geom_sf(data = POVERTY65.slope, aes(fill = score)) + 
  labs(title = "Poverty 65+ Score, New Jersey") +
  mapTheme() + 
  scale_fill_gradient2("Score", limits = c(1, 2), labels=c("1","2"),breaks=c(1,2),low = "#762A83", high = "#1B7837")

3.6 Combine All Scores

All scores are combined into one dataframe. The dataframe is then set to all values greater than “1” to “0”. This leaves a file where all grid cells that could benefit from affordable housing have a value of “1” and all others have a value of “0”. The data is then saved to a csv file for use in the SHINY application. Below is a facetted map series of the output.

knitr::opts_chunk$set(cache=TRUE)

#Combine all scores to one dataframe
study.panel <- fishnet
study.panel$AFFORDABILITY <- AFFORD.slope$score[match(study.panel$uniqueID, PROPVAL.slope$uniqueID)]
study.panel$PROPERTYVALUE <- PROPVAL.slope$score[match(study.panel$uniqueID, PROPVAL.slope$uniqueID)]
study.panel$DIVERSITY <- PCTWHITE.slope$score[match(study.panel$uniqueID, PCTWHITE.slope$uniqueID)]
study.panel$TRAVELTIME <- TRAVEL.slope$score[match(study.panel$uniqueID, TRAVEL.slope$uniqueID)]
study.panel$RETAXES <- RE.slope$score[match(study.panel$uniqueID, RE.slope$uniqueID)]
study.panel$ASTHMA <- ASTHMA.slope$score[match(study.panel$uniqueID, ASTHMA.slope$uniqueID)]
study.panel$EDUCATION <- HIGHSCHOOL.slope$score[match(study.panel$uniqueID, HIGHSCHOOL.slope$uniqueID)]
study.panel$POVERTY65 <- POVERTY65.slope$score[match(study.panel$uniqueID, POVERTY65.slope$uniqueID)]

#Keep only scores of 1

study.panel$AFFORDABILITY[study.panel$AFFORDABILITY > 1] <- 0
study.panel$PROPERTYVALUE[study.panel$PROPERTYVALUE > 1] <- 0
study.panel$DIVERSITY[study.panel$DIVERSITY > 1] <- 0
study.panel$TRAVELTIME[study.panel$TRAVELTIME > 1] <- 0
study.panel$RETAXES[study.panel$RETAXES > 1] <- 0
study.panel$ASTHMA[study.panel$ASTHMA > 1] <- 0
study.panel$EDUCATION[study.panel$EDUCATION > 1] <- 0
study.panel$POVERTY65[study.panel$POVERTY65 > 1] <- 0

#Save as CSV and export to SHINY Application Folder
study.panel.csv <- study.panel
study.panel.csv$geometry <- NULL
write.csv(study.panel.csv,"C:/Users/scott/Desktop/MUSA800/SHINY/appData.csv", row.names = FALSE)

#Facet Map
study.panel.Map <- study.panel %>%
  gather(Variable, Score, 4:11)

study.panel.Map$Score[study.panel.Map$Score == 0] <- NA

ggplot() + 
  geom_sf(data = study.panel.Map, aes(fill = Score)) + 
  labs(title = "Benefit from affordable housing") +
  mapTheme() + 
  facet_wrap(~Variable, ncol = 4)  + 
  scale_fill_gradientn(colours = "darkgreen")

#Export Fishnet to SHINY Application Folder
#setwd("C:/Users/scott/Desktop/MUSA800/SHINY") 
#st_write(fishnet, "fishnet.shp")
#setwd("C:/Users/scott/Desktop/MUSA800")

4.1 SHINY Application

+INSERT SHINY APPLICATION CODE+

4.2 SHINY Application URL

+INSERT SHINY APPLICATION URL+

4.3 SHINY Application Screenshot

+INSERT SHINY APPLICATION SCREENSHOT+

5.1 Conclusion

+INSERT CONCLUSION+